Fusariose: QTLs Description - Part2

This file follows the first part of the QTL analysis (Description_of_QTLs-Part1). Load the results + functions of the first part:

load("/Users/yan/Dropbox/Publi_Fusariose/ANALYSIS_REPRO/DATA/QTL_R_environment.R")

Charge some libraries that will be useful

library(RColorBrewer)
library(xtable)
library(plotly)
library(FactoMineR)

7/ IC on the map (Fig)

We can also visualize the IC on the genetic map. Let’s visualize important variable ONLY = AUDPC and last note for BLUPs only.

par(mfrow=c(1,1))
my_var_to_show=PUBLI
show_IC_QTL_on_map(data, my_var_to_show , c("1B", "5A"), my_QTL_thres , 1.5)

I save it as a png image for publication purpose.

png("~/Dropbox/Publi_Fusariose/FIGURE/Fig_IC_QTL_on_map.png")
show_IC_QTL_on_map(data, my_var_to_show , c("1B", "5A"), my_QTL_thres , 1.5)
dev.off()
## quartz_off_screen 
##                 2

8/ Evolution of LOD scores with time?

We are going to study qtl 1B and 5A only. Let’s have a look to the maximum LOD score of these QTL at each date of observation. For each experiment, for each variable?

Functions

Fonction qui regroupe chaque LOD max dans un tableau. On ne garde que les BLUPs de toutes les variables à toutes les dates.

find_LOD_max=function( chromo, my_min, my_max   ){
  bilan=data.frame(matrix(0, 200  ,8))
  num=0
  vec_names=c("NOT", "PEPI", "NEPIL", "PEPIL")
  for(expe in list(CAP13, GRI11, GRI13, GRI15, LEC14, BEQ11, BAR14)){
    
    num_carac=0
    for(carac in list(NOT, PEPI, NEPIL, PEPIL)){
      num_carac=num_carac+1
      my_carac=intersect(BLUP, intersect(expe, carac ))
      if( expe==BEQ11){ my_carac=intersect(expe, carac ) }
      if( length(my_carac)==0) { next }
      num=num+1
      my_vec=c()
      for(i in c(1:length(my_carac))){
        b=(data$LOD[ which(data$LG==chromo & data$variable==my_carac[i] & data$Distance<my_max & data$Distance>my_min ) ])
        a=max( b , na.rm=T)
        my_vec=c(my_vec, a)
      }
      
      my_replace=c( unique(gsub("_.*","",my_carac)) , vec_names[num_carac], my_vec )
      bilan[num, c(1:length(my_replace)) ] = my_replace 
    }
  }
  bilan=bilan[bilan[,1]!=0 , ]
  bilan[bilan==0]=NA
  return(bilan)
}

Fonction pour visualiser ce tableau.

show_max_LOD=function(bilan){
  
  # Préparation des axes des X.
  x1=c("300", "400", "550", "AUDPC")
  x2=c("300", "350", "450", "500", "550", "AUDPC")
  x3=c("350", "450", "550", "AUDPC")
  x4=c("350", "550", "AUDPC")
  x5=c("300", "400", "550", "AUDPC")
  x6=c("300", "350", "450", "500", "AUDPC")
  x7=c("300", "400", "550", "600", "AUDPC")
  my_x_axis=list(x1,x2,x3,x4,x5,x6,x7)
  
  par(mfrow=c(2,4) , mar=c(3,3,1,1))
  liste_var=c("NOT", "PEPI","NEPIL", "PEPIL")
  num=0
  for(expe in list("CAP13", "GRI11", "GRI13", "GRI15", "LEC14", "BEQ11", "BAR14")){
    num=num+1
    AA=bilan[ which(bilan[,1]==expe),]
    plot(seq(1,6) , AA[ 1, c(3:ncol(AA))] , type="l", ylim=c(0,10) , col=my_colors[match( AA[ 1, 2] , liste_var )], lwd=2, xlab=expe, ylab="", xaxt="n")
    axis(1, at=seq( 1 , length(my_x_axis[[num]])), my_x_axis[[num]], cex=0.7)
    abline(h=my_QTL_thres, col="grey")
    abline(v=c(1:6), col="grey", lwd=0.3)
    for(i in c(2:nrow(AA))){
     points(seq(1,6) , AA[ i, c(3:ncol(AA))] , type="l" , col=my_colors[match( AA[ i, 2] , liste_var )], lwd=2)  
    }
   text(x=1.8, y=9, expe, col="orange" )
  }
  plot(1,1,col="transparent", xaxt="n", yaxt="n", bty="n", xlab="", ylab="")
  legend( "center", legend = c("NOT", "PEPI","NEPIL", "PEPIL") , bty="n", pch=20 , col=my_colors[1:4], pt.cex=4, cex=1.6 )
}

QTL 1B

Calcul des LOD max à chaque date:

max_lod_1B=find_LOD_max( "1B", 0, 9 )

Show it:

show_max_LOD(max_lod_1B)

Donc ces graphiques représentent les valeurs des LODs max de chaque QTL, i.e. pour chaque variable de chaque expe pour le chromosome 1B.

QTL 5A

Calcul des LOD max à chaque date:

max_lod_5A=find_LOD_max( "5A", 100, 300   )

Show it:

show_max_LOD(max_lod_5A)

10/ Comparaison avec les QTLs de la biblio

Using SSR

Let’s have a look to all the 9 SSR markers I have in my map:
Is it normal I have only 9?
This is not gonnna help a lot…
- 1B: pas de SSR
- 5A: 2 SSR, mais très très loin de l’IC.

–> Il faut aller voir si je peux trouver des arqueurs de la littérature qui ont été passé sur Dic2 x Silur, et les replacer sur la carte comme pour la publi mosaiquE…

SSR=unique(data[ grep("gwm" , data$marqueur) , c(1:3) ])
nrow(SSR)
## [1] 9
LG marqueur Distance
1A Xgwm164 68.40
3B Xgwm376 92.30
3B Xgwm77 94.20
4B Xgwm495 94.00
5A Xgwm205 41.00
5A Xgwm415 72.50
5B Xgwm544 55.30
7A Xgwm276 181.50
7B Xgwm400 38.50

Position of Snn1 on map

Le gène Snn1 est positionné sur le 1B et apporte une résistance forte à la septoriose, un autre champignon. On peut trouver la séquence de ce gène ici:
https://www.ncbi.nlm.nih.gov/nuccore/858936371/ Publié par: http://advances.sciencemag.org/content/2/10/e1600822.full

Blast sur la référence blé tendre release 28

makeblastdb -in sequence_snn1.fasta -dbtype nucl -out database

Blast à 100% sur ce contig: Traes_1BS_C4B727779

Ce contig n’est pas dans notre carte génétique, donc on a pas de SNP dessus. Par contre, en terme de position physique, il se situe à 880274 Mb sur le 1BS. Ca correspondrait donc à être entre 1.19 et 1.59 cM sur notre carte.

1B Traes_1BS_7BC6EFB55@1326 1.19999999999999 1B 502322 1B Traes_1BS_81F7DAE3E@432 1.59999999999999 1B 2538668

Nous notre QTL a tendance a etre un poil plus loin, mais vraiment vraiment proche…

11/ Relationship QTL presence / field characteristic

The intensity of the Fusarium infection was not the same from one experiment to another. Let’s check the relationship between disease intensity and strength of the QTL detected?

12/ Epistatic effect?

OK so we have 2 main QTLs. Do we have a epistatic interaction betwwen them??

Example

Let’s start on an example: BAR14_PEPIL300_blup. 2 QTLs: markers Cluster_16778|Contig1|original@1840 et Cluster_9940|Contig1|complementarySeq@104

# calculate inter
bil_inter=data.frame(matrix(0,0,9))
colnames(bil_inter)=c("trait","marker-1B", "marker 5A", "mean A-A", "mean A-B", "mean B-A", "mean B-B", "R2 tot", "pval inter")
bil_inter[1,]=analyse_inter("Cluster_16778|Contig1|original@1840", "Cluster_9940|Contig1|complementarySeq@104" , "BAR14_PEPIL300_blup")
trait marker-1B marker 5A mean A-A mean A-B mean B-A mean B-B R2 tot pval inter
BAR14_PEPIL300_blup Cluster_16778|Contig1|original@1840 Cluster_9940|Contig1|complementarySeq@104 -0.2 3.32 -0.91 -0.65 0.281 0.0001369779

And show this interaction!

boxplot_two_QTL_interactive("Cluster_16778|Contig1|original@1840", "Cluster_9940|Contig1|complementarySeq@104" , "BAR14_PEPIL300_blup")
#boxplot_two_QTL("Cluster_16778|Contig1|original@1840", "Cluster_9940|Contig1|complementarySeq@104" , "BAR14_PEPIL300_blup", ylim=c(-4,7))

All trait

Il faut commencer par trouver tous les traits qui ont un QTL signif sur le 1B et un signif sur le 5A.

var_with_2QTLs=c()
for(i in list(summary_QTL_NOT, summary_QTL_PEPI, summary_QTL_NEPIL, summary_QTL_PEPIL)){
  a=i[which(i$chromo%in%c("1B","5A")) , ]
  b=table(a$carac)
  var=names(b[b==2])
  var_with_2QTLs=c(var_with_2QTLs, var)
}
var_with_2QTLs=var_with_2QTLs[grep("blup",var_with_2QTLs)]
On a 8 variables qui ont les 2 QTLs significatifs. Voila de quelles variables il s’agit:
var_with_2QTLs
BAR14_PEPI350_blup
BAR14_PEPI400_blup
BAR14_PEPI600_blup
BAR14_PEPIAUDPC_blup
GRI13_PEPI550_blup
BAR14_NEPIL300_blup
BAR14_PEPIL300_blup

Calculons les p-values de l’effet d’interaction pour tous ces cas:

for(i in var_with_2QTLs[-which(var_with_2QTLs=="BAR14_PEPIL300_blup")]){
  vv=analyse_inter("Cluster_16778|Contig1|original@1840", "Cluster_9940|Contig1|complementarySeq@104" , i)
  bil_inter=rbind(bil_inter , vv)
}
trait marker-1B marker 5A mean A-A mean A-B mean B-A mean B-B R2 tot pval inter
BAR14_PEPIL300_blup Cluster_16778|Contig1|original@1840 Cluster_9940|Contig1|complementarySeq@104 -0.2 3.32 -0.91 -0.65 0.281 0.0001369779
BAR14_PEPI350_blup Cluster_16778|Contig1|original@1840 Cluster_9940|Contig1|complementarySeq@104 0.38 21.58 -7.91 -3.92 0.355 0.0004079017
BAR14_PEPI400_blup Cluster_16778|Contig1|original@1840 Cluster_9940|Contig1|complementarySeq@104 -0.08 20.57 -7.65 -4.53 0.244 0.0039045779
BAR14_PEPI600_blup Cluster_16778|Contig1|original@1840 Cluster_9940|Contig1|complementarySeq@104 2.17 3.69 -4.75 1.02 0.318 0.0163055323
BAR14_PEPIAUDPC_blup Cluster_16778|Contig1|original@1840 Cluster_9940|Contig1|complementarySeq@104 896.84 9825.26 -4380.44 -2080.69 0.394 0.002065465
GRI13_PEPI550_blup Cluster_16778|Contig1|original@1840 Cluster_9940|Contig1|complementarySeq@104 -1.25 0.94 0.54 1.12 0.023 0.4583633636
BAR14_NEPIL300_blup Cluster_16778|Contig1|original@1840 Cluster_9940|Contig1|complementarySeq@104 -1.33 18.18 -4.85 -3.36 0.259 0.0002526347

On a donc un effet d’épistasie qui est significatif à tous les coups (quasiment).

Représentons ces épistasies avec un boxplot pour chacune des variables:

# widget to choose data
#selectInput("select", label = h3("Select a variable"), choices = var_with_2QTLs,  selected = 1)

# interactibe boxplot
#renderPlotly({
#  boxplot_two_QTL("Cluster_16778|Contig1|original@1840", "Cluster_9940|Contig1|complementarySeq@104" , select)
#})

Dans la majorité des cas, on observe que les individus avec 1 seul QTL ou avec 2 QTLs sont similaires. Par contre, les individus avec aucun des 2 sont très différents et très sensibles. Donc en gros, si j’ai au moins un des 2 gènes tout va bien, sinon je suis sensible. Donc on est pas dans un modèle additif à priori, les 2 gènes sont liés.

Yan Holtz

October 2016